{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Check if all the packages are installed or not\n",
    "cond = \"Gadfly\" in keys(Pkg.installed()) &&\n",
    "\"Colors\" in keys(Pkg.installed()) &&\n",
    "\"ODEInterface\" in keys(Pkg.installed()) &&\n",
    "\"ForwardDiff\" in keys(Pkg.installed());\n",
    "@assert cond \"Please check if the following package(s) are installed:\\\n",
    "    Gadfly\\\n",
    "    Colors\\\n",
    "    ODEInterface\\\n",
    "    ForwardDiff\"\n",
    "\n",
    "# Load all the required packages\n",
    "using Gadfly\n",
    "using Colors\n",
    "using ODEInterface\n",
    "using ForwardDiff\n",
    "@ODEInterface.import_huge\n",
    "loadODESolvers();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "######################## Function for saving plots #################################\n",
    "# Input:\n",
    "# fileName = Name of the file where the plot is to be stored\n",
    "#            (with or without extension)\n",
    "# f_e = Array containing function evaluations as columns for each solver\n",
    "# err = Array containing erros as columns for each solver\n",
    "# solverNames = Array containing the names of solvers used in respective order\n",
    "# plotSize = size of the plot to be created\n",
    "###################################################################################\n",
    "function savePlotPNG(fileName,f_e,err,solverNames,\n",
    "    plotSize=[30cm,30cm])\n",
    "    \n",
    "    numOfLayers = length(solverNames);\n",
    "    \n",
    "    if !contains(fileName,\".\")\n",
    "        fileName = string(fileName,\".png\");\n",
    "    end\n",
    "    \n",
    "    plotColorsHex = [\"#4D4D4D\",\"#5DA5DA\",\"#FAA43A\",\"#60BD68\",\n",
    "        \"#F17CB0\",\"#B2912F\",\"#B276B2\", \"#DECF3F\",\"#F15854\"];\n",
    "    plotColors = [parse(Colorant,c) for c in plotColorsHex];\n",
    "    \n",
    "    majorFontSize = 24pt;\n",
    "    minorFontSize = 20pt;\n",
    "    pointSize = 5pt;\n",
    "    \n",
    "    myplot = plot(Scale.x_log10,Scale.y_log10,\n",
    "        Coord.cartesian(xflip=true,xmin=-11),\n",
    "        Guide.manual_color_key(\"Legend\",solverNames,plotColorsHex[1:numOfLayers]),\n",
    "    Guide.xlabel(\"error\"),Guide.ylabel(\"#Function evaluations\"),\n",
    "        Theme(major_label_font_size=majorFontSize,panel_stroke=colorant\"black\",\n",
    "        minor_label_font_size=minorFontSize,key_title_font_size=majorFontSize,\n",
    "        key_label_font_size=minorFontSize,key_position=:top,key_max_columns=1));\n",
    "    \n",
    "    for i = 1:numOfLayers\n",
    "        push!(myplot,layer(x=err[:,i],y=f_e[:,i],Geom.point,Geom.path,\n",
    "        Theme(default_color=plotColors[i],default_point_size=pointSize)));\n",
    "    end\n",
    "    \n",
    "    draw(PNG(fileName,plotSize[1],plotSize[2]),myplot)\n",
    "    return nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define the right-hand function for automatic differentiation\n",
    "function vdpolAD(x)\n",
    "    return [x[2],((1-x[1]^2)*x[2]-x[1])*1e6]\n",
    "end\n",
    "\n",
    "# Define the system for the solver\n",
    "function vdpol(t,x,dx)\n",
    "    dx[:] = vdpolAD(x);\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "# Define the Jacobian function using AD\n",
    "function getJacobian(t,x,J)\n",
    "    J[:,:] = ForwardDiff.jacobian(vdpolAD,x);\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "# Flag to check whether plot is to be generated and saved or not\n",
    "# Also checks if all solvers are successful\n",
    "printFlag = true;\n",
    "\n",
    "# Initial conditions\n",
    "t0 = 0.0; T = [1.0:11.0;]; x0 = [2.0,0.0];\n",
    "\n",
    "# Get \"reference solution\" from\n",
    "# http://www.unige.ch/~hairer/testset/testset.html\n",
    "f = open(\"vdpolRefSol.txt\");\n",
    "lines = readlines(f);\n",
    "\n",
    "x_ref = Array{Float64}(11);\n",
    "\n",
    "tmp = Array{Float64}(22);\n",
    "counter = 1;\n",
    "for l in lines\n",
    "    tmp[counter] = parse(Float64,l);\n",
    "    counter += 1;\n",
    "end\n",
    "\n",
    "x_ref = tmp[1:2:end];\n",
    "\n",
    "close(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Store the solver names for plotting\n",
    "solverNames = [\"RADAU\",\"RADAU5\",\"SEULEX\"];\n",
    "\n",
    "# Initialize the variables for plots\n",
    "# f_e = number of function evaluations\n",
    "f_e = zeros(33,3);\n",
    "# err = error wrt ref solution over all time steps and components  \n",
    "err = zeros(33,3);\n",
    "# flops = Floating point operations\n",
    "flops = zeros(33,3);\n",
    "\n",
    "# Weights for computing flops\n",
    "dim = 2; # dimension of the system\n",
    "flopsRHS = 5; # Counted\n",
    "flopsLU = ceil(2*((dim)^3)/3); # As per LU algorithm (can be less)\n",
    "flopsFW_BW = (dim)^2; # As per FW/BW algorithm (can be less)\n",
    "flopsJac = ceil(1.5*flopsRHS); # A guess at the moment\n",
    "\n",
    "for i =0:32\n",
    "\n",
    "    # Set the tolerance for current run\n",
    "    Tol = 10^(-2-i/4);\n",
    "\n",
    "    # Set solver options\n",
    "    opt = OptionsODE(OPT_EPS=>1.11e-16,OPT_ATOL=>Tol,OPT_RTOL=>Tol,\n",
    "    OPT_RHS_CALLMODE => RHS_CALL_INSITU,\n",
    "    OPT_JACOBIMATRIX=>getJacobian);\n",
    "\n",
    "    # Store the stats of the last t_end\n",
    "    # for computing flops\n",
    "    stats = Dict{ASCIIString,Any};\n",
    "\n",
    "    # Restart the solution for each end time\n",
    "    # to ensure a more accurate solution\n",
    "    # compared to dense output\n",
    "\n",
    "    # Solve using RADAU\n",
    "    x_radau = Array{Float64}(11);\n",
    "    for j=1:11\n",
    "        (t,x,retcode,stats) = radau(vdpol,t0, T[j], x0, opt);\n",
    "        # If solver fails do not continue further\n",
    "        if retcode != 1\n",
    "            println(\"Solver RADAU failed\");\n",
    "            printFlag = false;\n",
    "            break;\n",
    "        end\n",
    "        x_radau[j] = x[1];\n",
    "        f_e[i+1,1] = stats.vals[13];\n",
    "    end\n",
    "    # If solver fails do not continue further\n",
    "    if !printFlag\n",
    "        break;\n",
    "    end\n",
    "    err[i+1,1] = norm(x_radau-x_ref,Inf);\n",
    "    flops[i+1,1] = stats[\"no_rhs_calls\"]*flopsRHS+\n",
    "                   stats[\"no_fw_bw_subst\"]*flopsFW_BW +\n",
    "                   stats[\"no_jac_calls\"]*flopsJac+\n",
    "                   stats[\"no_lu_decomp\"]*flopsLU;\n",
    "\n",
    "    # Solve using RADAU5\n",
    "    x_radau5 = Array{Float64}(11);\n",
    "    for j=1:11\n",
    "        (t,x,retcode,stats) = radau5(vdpol,t0, T[j], x0, opt);\n",
    "        # If solver fails do not continue further\n",
    "        if retcode != 1\n",
    "        println(\"Solver RADAU5 failed\");\n",
    "            printFlag = false;\n",
    "            break;\n",
    "        end\n",
    "        x_radau5[j] = x[1];\n",
    "        f_e[i+1,2] = stats.vals[13];\n",
    "    end\n",
    "    # If solver fails do not continue further\n",
    "    if !printFlag\n",
    "        break;\n",
    "    end\n",
    "    err[i+1,2] = norm(x_radau5-x_ref,Inf);\n",
    "    flops[i+1,2] = stats[\"no_rhs_calls\"]*flopsRHS+\n",
    "                   stats[\"no_fw_bw_subst\"]*flopsFW_BW +\n",
    "                   stats[\"no_jac_calls\"]*flopsJac+\n",
    "                   stats[\"no_lu_decomp\"]*flopsLU;\n",
    "\n",
    "    # Solve using SEULEX\n",
    "    x_seulex = Array{Float64}(11);\n",
    "    for j=1:11\n",
    "        (t,x,retcode,stats) = seulex(vdpol,t0, T[j], x0, opt);\n",
    "        # If solver fails do not continue further\n",
    "        if retcode != 1\n",
    "            println(\"Solver seulex failed\");\n",
    "            printFlag = false;\n",
    "            break;\n",
    "        end\n",
    "        x_seulex[j] = x[1];\n",
    "        f_e[i+1,3] = stats.vals[13];\n",
    "    end\n",
    "    # If solver fails do not continue further\n",
    "    if !printFlag\n",
    "        break;\n",
    "    end\n",
    "    # Get the error over all the components and\n",
    "    err[i+1,3] = norm(x_seulex-x_ref,Inf);\n",
    "    flops[i+1,3] = stats[\"no_rhs_calls\"]*flopsRHS+\n",
    "                   stats[\"no_fw_bw_subst\"]*flopsFW_BW +\n",
    "                   stats[\"no_jac_calls\"]*flopsJac+\n",
    "                   stats[\"no_lu_decomp\"]*flopsLU;\n",
    "end\n",
    "\n",
    "# Save the plot in PNG format\n",
    "# if all the solvers were successful\n",
    "if printFlag\n",
    "    savePlotPNG(\"vdpolPrecisionTest\",f_e,err,solverNames);\n",
    "else\n",
    "    println(\"Plot cannot be generated due to failure\");\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Old reference solution\n",
    "Tol = 1e-14;\n",
    "# for Tol < 1e-14 we get the error \"TOLERANCES ARE TOO SMALL\"\n",
    "opt = OptionsODE(OPT_EPS=>1.11e-16,OPT_RTOL=>Tol, OPT_ATOL=>Tol,\n",
    "OPT_RHS_CALLMODE => RHS_CALL_INSITU,\n",
    "OPT_JACOBIMATRIX => getJacobian);\n",
    "\n",
    "# Store only the desired component\n",
    "# Here, only the first component is desired\n",
    "# The second component is the first derivative of the first component\n",
    "# due to the fact that it is a second order system.\n",
    "# Hence error will be taken over the first component only.\n",
    "x_ref = Array{Float64}(11);\n",
    "\n",
    "for i=1:11\n",
    "    (t,x,retcode,stats) = seulex(vdpol,t0, T[i], x0, opt);\n",
    "    if retcode!=1\n",
    "        printFlag = false;\n",
    "        break;\n",
    "    end\n",
    "    x_ref[i] = x[1];\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
